Determinations of upper critical field in continuous Ginzburg-Landau model 
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Novel procedures to determine the upper critical field Bc2 have been proposed within a continuous 
Ginzburg-Landau model. Unlike conventional methods, where Bc2 is obtained through the determi- 
nation of the smallest eigenvalue of an appropriate eigen equation, the square of the magnetic field 
is treated as eigenvalue problems so that the upper critical field can be directly deduced. The calcu- 
lated Bc2 from the two procedures are consistent with each other and in reasonably good agreement 
with existing theories and experiments. The profile of the order parameter associated with Bc2 is 
found to be Gaussian-like, further validating the methodology proposed. The convergences of the 
two procedures are also studied. 



I. INTRODUCTION 

The determinations of the upper critical field Bc2 are 
based on the Ginzburg-Landau- Abrikosov-Gorkov-'^-^-'^ 
(GLAG) framework that is applicable to practically all 
superconductors. The starting point of the macro- 
scopic description is the first Ginzburg-Landau (GL) 
equationi*^ or its similarities. Examples would be 
the standard linearized GL equation^, the harmonic 
oscillator equationSiSiSii, the Mathieu equatior&SiiS, 
the continuous Ginzburg-Landau equationii, the GLd 
equation(s)^^-^''-^'^-^^ [Ginzburg-Landau model for (d-|-s)- 
wave superconductors], and the Ginzburg-Pitaevskii 
equationi^. 



On the other hand, the linear gap equatio] 
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adopted to microscopically determine the upper critical 
field. This equation can be converted into an equation 
including the digamma function. The parameter of the 
digamma function is related to the smallest eigenvalue of 
an eigen equation with appropriate boundary conditions. 
Having obtained the smallest eigenvalue of the eigen 
equation, Bc2 may then be implicitly determined from 
the equation including the digamma function {Bc2 is 
in connection with the smallest eigenvalue). Such a 
procedure or its extensions may be seen, for example, in 
the calculations of i?c2 for type-II superconductorsiSiSS 
(s-wave), heavy-fermion superconductorsSiS^ (p- 
wave), high Tc superconductora^2iS4 [(d-l-s)-wave] , 
layered superconductor a^°i^^i^^ , superlattice a^"^!^^ , thin 
fllms22i2£, multilayers^ and organic superconductors^. 



It is tempting to conclude that almost all the 
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made to obtain Bc2 are via the determination of the 
smallest eigenvalue of an appropriate eigen equation. 
For instance, in Ref. 0, Bc2 is inversely proportional to 
the smallest eigenvalue (the lowest Landau level) of the 
harmonic oscillator equation. 

Nevertheless, for all the GL-like equations as well as 
the linear gap equation, one may have a universal defini- 
tion for the upper critical field: it is the maximum mag- 
netic field at which the corresponding equation has a non- 
trivial solution^. Such a nontrivial solution would nat- 
urally correspond to the eigen function associated with 



the smallest eigenvalue and would correspond to the so- 
lution to the linear gap equation at Bc2^'^- It is worth 
noting that the eigen function associated with the small- 
est eigenvalue may be just one of the solutions to the 
linear gap equation (see, for example, Refs. 0, and 

ii. 

Several approaches, for example, the perturbation 
methodSiiSii^ and the variational tcchnique^'^'^'*, have 
been adopted to find Bc2 by obtaining the smallest eigen- 
value. We note that in Ref. 4, the authors iteratively 
calculated the upper critical field of a thin-film supercon- 
ductor with a ferromagnetic dot: first guessing a value 
for the upper critical field, then solving the standard lin- 
earized GL equation to make sure the calculated value of 
the upper critical field is equal to the guessed one. Oth- 
erwise, equate the newly calculated value to the guessed 
one and resolve the linear equation until the two values 
are equal to each other. It can be seen that such an it- 
erative method is still related to obtaining the smallest 
eigenvalue"* . 

In this paper, instead of obtaining the smallest eigen- 
value, from which Bc2 is implicitly or indirectly deter- 
mined, we will directly calculate Bc2- the square of the 
magnetic field B^ is treated as an eigenvalue of an eigen 
equation so that Bc2 can be directly deduced from the 
largest eigenvalue of the eigen equation. Within a contin- 
uous Ginzburg-Landau model^, two procedures will be 
presented to obtain the corresponding eigen equations, 
either of which can determine Bc2- Note that the applica- 
tions of these procedures have been published in Refs. 3^ 
andlM 



II. MODEL 

The continuous Ginzburg-Landau model (CGL) is ap- 
plicable to layered superconductors^^. The unit cell de- 
scribing the layering structure consists of alternating su- 
perconducting (S) and insulating (I) layers. The z axis 
is normal to the layers and its origin is at the midpoint 
of one of the I layers. The center of the S layer is located 
at D/2, where D is the size of the unit cell^. The GL 
coefficients and the effective superpair masses (perpen- 
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dicular and parallel to the layers) in the CGL free energy are spatially dependent. The CGL free energy is, 
I 



F 



df / dz 



a{T, z)|vl/(r, z)\^ + l/3|*(f, + 
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2m{z) 



_'^A^-'\r,z)]^{r,z) 



l^B^r^z) 



(1) 



where T is the temperature, r = (x, y) is the planar 
vector and A{f,z) = {A'^'^\r, z), Az{f, z)) is the vector 
potential. The effective masses and the GL condensa- 
tion coefficient a{T, z) are assumed to be periodic with 
a period D. j3, however, is held fixed as it does not af- 
fect the qualitative behavior of the studied system'^^ . As 
before^SiSLSS^a, we choose 



a(T, z) 
1 



M{z) 
1 

m(z) 



[ao + aiCOs{2TTz/D)]{l-T/T,), (2a) 
Go-}-GiCos(27rz/L>), (2b) 



.90 + .91 cos(27rz/£'). 



(2c) 



where Tc is the transition temperature, ao, ai, Gq, Gi, 
go and gi are model parameters. Note that for dirty 
materials, the linear GL theory for Bc2 may be extended 
to low temperatures^iiSiS^. 

Let an external field B be applied along the y-direction 
so that the vector potential can be taken as A = 
iBz,0,0). Assuming ^'(f, z) = e''=ii-^*(z) and minimiz- 
ing the free energy of Eq. ^ we get 
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(3) 



with Zs = hkx/{2eB). At B = Bc2, the order param- 
eter is small enough so that the term /3|^(0)p\I'(2) in 
Eq. O can be omitted. The superconducting order at 
Bc2 nucleates in the S layer first so that one may choose 
Zs — I?/3ii. To explore the features of the order pa- 
rameter along the z-direction, one may assume ky = 0. 
Finally, we arrive at the following equation. 
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At a given temperature T, the maximum magnetic field B 
at which a nontrivial solution satisfies the above equation 
gives a point in the Bc2-T plot. Eq.^will be numerically 
solved subject to the following boundary conditions. 



*(0) = *(i5). 



d_ 

dz 



(5a) 
(5b) 



z=D 



Note that neglecting the size effect in the z-direction and 
considering the layering structure and the superposition 
of different superpairs as a mean field effect, the variation 
of the macroscopic wave function-the order parameter 
along the z-direction may have a periodic property (for 
example, see Ref. IlL) . It should be mentioned that the 
new methodology of determining the upper critical field, 
which is to be described below, is also applicable to other 
boundary conditionsSSiiS. 



III. METHODOLOGY 



Procedure I 



Taking into account the boundary conditions of Eq. |31 
Eq. 21 can be transformed into a system of equations 
which can be simplified to 



U* = 0, 



(6) 



where ^ = ('Si, ^2, ^2n+i)' is a column vector rep- 
resenting the discrete solutions of Eq. 0] The symbol ' 
here indicates transpose and 71 is a positive integer. U is 
a {2n + 1) X {2n + 1) sparse matrix having the following 
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where Ai^ Bt, d, Di and Ei are the coefhcients of the 
discretized equations of Eq. ^ In the first row of U, the 
boundary condition ^'(0) = ^{D) CEa. I5b|) is exphcitly 
written as ^'i = 4'2n+i while in the last row, the dis- 
cretization of -§;'^{z)\^^u — (Eg. I5b|l is implemented 
as -2«'2 - *2«-i + 6*^„ - 3^'2„+i = 0. The periodic 
property of the solutions is taken into consideration in 
the last row (5'2n+2 = ^2), in the second row (see A2) 
and in the second last row (see Em). 

It is easy to verify that only the Ci coefficients con- 
tain the magnetic field B and these coefficients can 
be separated into two terms, where one of them is B- 
dependenlj^, 

To obtain B^i with our method, we require that there 
be no prefactor appearing before the square of the mag- 
netic field B^ . Hence, the factor Cf in Eq. |H1 should 
be divided. To cure the division singularity at z = -D/2 
(i.e., cf = Oati = n+ lin Eq. O, a solution property 
^5'(z)|^_^,2 = is assumed and implemented in the 



middle row of U as -2\E'„ - 3*„+i -f 6\I'„+2 - *,i+3 = 
(different from the approximation in the last row). 

For Eq.|Hlto have non-trivial solutions, the determinant 
of U should be zero, 

det|U| = 0. (9) 

By eliminating the constant elements in the first, middle 
and last rows of U and separating Ci into the two terms 
with the negative prefactor —Cf excluded (cf. Eq. |HJ|, 
det|U| can be transformed into 

det|P - B^ii = 0, (10) 
where I is a unitary matrix and P is a matrix without the 
magnetic field (see below). Thus, the largest solution for 
B, namely i?c2, can be easily available just by obtaining 
the largest positive eigenvalue of the eigen equation, 

P = B2$, (11) 

where $ is the eigen function and P, with the dimension 
of (2n — 2) X (2n — 2), has the following structure. 
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where A', = ~A,/Cf , i?^ = -B./Cf, Cf - -C^o/Cf , 
= -D,/CP and i;,' = -E./Cf. Here, i 7^ n + 1 and 
thus the division singularity is avoided. Having deter- 
mined i?c2 from Eg. 111! one can obtain the corresponding 
order parameter by substituting 3^2 back into Eq. 1^1 

B. Procedure II 

In the above procedure, the upper critical field is ob- 
tained by treating the square of the magnetic field 



as an eigenvalue (Eq. . Following the same idea, one 
may obtain a straight forward procedure from Eq. ^ 

Q* = S2*, (13) 

where 'S' — (^'i, 5*2, ^2n)', which has the same mean- 
ing as that in Eq. However, to avoid the singularity 
at z = D/2, the dimension of the current 4*, 2n, is set 
different from the previous one, 2n + 1. The matrix Q 
has the following structure^, 
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(14) 



The dimension of Q is 2n x 2n, larger than that of P but Eq. 1131 implies that one may directly determine the 
smaller than that of U. upper critical field from the linear GL equation. The 
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present procedure is essentially consistent with the usual 
method for obtaining B(.2 through the smallest eigenvalue 
of the linear GL equation since there exists an inverse 
relation between the magnetic field eigenvalues and the 
eigenvalues (or the indices of the Landau levels) of the 
linear GL equation (for example, see Ref. 0). Our proce- 
dures (particularly Procedure II) indicate that one may 
directly determine the upper critical field without con- 
sidering the smallest eigenvalue of an appropriate eigen 
equation. In fact, in one of the earliest papers^, Gor'kov 
simply used a variational method to obtain 3^2- 



IV. RESULTS AND DISCUSSION 

Bi2212 has been chosen as our modeling prototype 
since it possesses a large anisotropy**^, and is thus suit- 
able for studies of properties related to the layered struc- 
ture. The values of the parameters for Bi2212 are taken 
from Ref . [sj The calculated i?c2 at zero temperature for 
various n are shown in Tabled It can be seen that the 
two sets of data calculated from Eas. II II andll3l converge 
finally and the convergent values are in good agreement 
with each other. 

A glance at the data for 5^2 in Table shows that 
Procedure II converges more rapidly than Procedure I. 
Hence, using Bc2{n = 1200) from Procedure II for the 
convergence analysis, the error is defined as 



En = \Bc2{n)- Bc2{n = l2m)\. 



(15) 



This equation is applicable to both procedures and the 
corresponding errors are listed in Tableland plotted in 
Fig. n It is found that the errors of the two procedures 
can be respectively fitted as 



En 
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for Procedure II. 



(16) 



TABLE L Values of Bc2 (SI unit) for various n, treated as 
eigenvalue problems from different procedures. The error is 
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FIG. 1: Errors for Procedure I and Procedure II, which are of 
order 3 and 4, respectively. 



Here no is taken to be 50 and En^ is the corresponding 
error. The grid spacing h ~ D/2n, when n is large. 
From these fits, we know that the error of Procedure I 
is of order 3 while that of Procedure II is 4, which are 
consistent with the error orders of the approximations 
used in the matrices U and Q, respectively. 

Note that Bc2 at K in the tables may be compara- 
ble to those extrapolated (using the WHH theory^") from 
some experiments (for example, see Ref. 43). By choosing 
appropriate values of the model parameters, the extrap- 
olated result of 2640 Tesla for Bi2212'*3 ^-.g obtained 
exactly. 

We find that the order parameters obtained from Eqs.El 
and 1 131 are also consistent with each other. In the follow- 
ing calculations, Procedure II is employed. In Fig. |21 we 
plot the spatial distribution of the calculated order pa- 
rameter at T/Tc = 0.9. It is found that for Bi2212 in a 
large temperature range, the asymptotic behavior of the 
order parameter can be expressed as 



^{z) ~ Cexp 



2e 



(17) 



which is exactly the ground state of the linearized GL 
equation at Bc2^-^ . The Gaussian profile of the order pa- 
rameter is also similar to the shape of the pair amplitude 
in bulk superconductors. In fact, it is natural that the 
order parameter of the largest eigenvalue of B^ should be 
just that of the lowest eigenvalue for the linear equation, 
due to the inverse relation between B and the eigenvalues 
of the linear equation^. 

When near Tc or in a less anisotropic superconductor 
such as YBCO, we found that an offset of the Gaussian 
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FIG. 2: Spatial distribution of the order parameter at T/Tc = 0.9, 
which can be approximated by the well-known Gaussian function. 



order parameter (Eq. 117(1 can not be neglected. Hence, 
one may infer that the mean field effect or the interlayer 
coupling'^® of the order parameter has less influence on 
a highly anisotropic superconductor at low temperatures 
(2-D feature) than on a superconductor at high temper- 
atures or with a small anisotropji^i. 

The finding here that a Gaussian-like exponential com- 
ponent contributes to the order parameter may ade- 
quately signify that the eigen function associated with 
the maximum magnetic field eigenvalue of the linear GL 
equation is equivalent, to a certain extent, to the eigen 
function related to the lowest eigenvalue of the linear GL 
equation. The latter eigen function may have the form 



exp 



(z - (cf. Eq. 4.73 in Ref. H). Compar- 

ing this function to the Gaussian-like exponential contri- 
bution (Eq.[T7|l, one has Bc2 oc 1/^^ p^j. Bi2212, we find 
that 1/^^ (X 1 — T/Tc is satisfied for a large temperature 



calculation 
• linear fit 




FIG. 3: Linear relation between 1/^^ and I — T/Tc, where § is the 
parameter characterizing the width of the Gaussian order parame- 
ter. 



FIG. 4: The Bc2-T behavior is linear in a large temperature range 
but square-root-like near Tc (inset). The solid circles (•) stand for 
the calculated data and the full lines are guidances to the eye. The 
dashed line represents a square-root fit. The reduced temperature 
of the transition point from square-root to linearity is about 0.997 
(indicated by the arrow). 



range (Fig.EJ, hence, 5^2 « 1 - T/T^ (see Fig.©. Sim- 
ilarly, we also have found that the linear Bc2-T relation 
is satisfied in YBCO. Such a linear feature is consistent 
with the general Bc2-T picture in layered cuprates near 
T = KM, and is in agreement with the anisotropic 
GL theory and some other theories2ii2iiS4&. Experimen- 
tally, the linear behavior away from Tc was observed in 

A closer inspection of Fig. ^reveals that near Tc, the 
Bc2-T plot exhibits a square-root-like behavior which is 
in agreement with Feinberg's theory^. This behavior 
is also found in twined superconductors*^, multilayers*-, 
and thin film a^i*^ . In one of the earliest papers, Ginzburg 
and Landaui correctly predicted that Bc2 for thin films 
varied as \f%^^^ when T — > Tc, as a consequence of the 
boundary conditions. This was subsequently confirmed 



experimentally by Khukhareva^". 

In our calculations, we found that near Tc, the height 
of the Gauss component C in Eq. El decreases with 
temperature and hence, the order parameter is nearly 
constant (3-D feature) ^^i*° . Qualitatively speaking, the 
nearly constant order parameter near Tc for the layered 
superconductors is equivalent to the case of thin films so 
that one can observe the square-root behavior near Tc- 
In fact, the boundary condition of Eq. I5bl does apply to 
thin films. More precisely, with a constant solution to the 
order parameter, one can immediately infer from Eq. 0| 
that there exists a square-root relationship between Bc2 
and 1 — T/Tc- It is worth noting that Dediu et al. also 
obtained a square-root Bc2 — T relation near Tc by solv- 
ing the usual GL equations and the feature of the order 
parameter near Tc was used to interpret their results. 

It is interesting to note that in the Bc2-T plot, the 
transition from square-root to linearity may be consid- 
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ered as a 3-D to 2-D crossover. The transition tempera- 
ture obtained here (indicated by the arrow in Fig. ^ is 
reasonably consistent with the estimates in Refs. ^ and 

m 

V. CONCLUSION 

The conventional work to determine upper critical 
fields should resort to obtaining the smallest eigenvalue 
of an appropriate eigen equation. However, within a con- 
tinuous Ginzburg-Landau model, we have demonstrated 
through two procedures that one can obtain the upper 
critical field by treating the square of the magnetic field 
as eigenvalue problems, from which Bc2 can be directly 
deduced. 

The calculated Bc2 from the two procedures are consis- 



tent with each other and in reasonably good agreement 
with existing theories and experiments. The profile of the 
order parameter obtained at Bc2 is Gaussian-like, further 
indicating the plausibility of the procedures proposed. 

The convergences of the proposed procedures were also 
investigated to identify the efficiency of the procedures. 
It was found that the more direct method, Procedure 
II, converges faster than Procedure I. The corresponding 
convergent error orders of the two procedures are 4 and 3, 
respectively. These orders are consistent with the orders 
of the approximations used in our calculations. 

Note that certain physical phenomena such as fluctua- 
tions and spin-orbit scattering may have some infiuences 
on the upper critical field. The present procedures pro- 
posed may serve as a starting point to further study more 
properties of upper critical fields. 
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